"""
Multivariate function model fitting with constraints
"""

import pandas as pd
import numpy as np
from scipy.optimize import minimize
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from matplotlib import pyplot as plt
import statsmodels.api as sm

plt.rcParams['font.sans-serif'] = ['Times New Roman']
plt.rcParams['axes.unicode_minus'] = False

data = pd.read_excel("归一化后的数据.xlsx")
data = data.fillna(0)

y = data.iloc[1:13, 4].values  # y is the fifth column, starting from the second row
x = data.iloc[1:13, :4].values  # x is from the second row, first to fourth columns
x1 = data.iloc[1:, :4].values  # x1 is from the second row onward, first to fourth columns

# Define the objective function for multivariate linear regression (least squares)
def objective(params, x, y):
    a = params[:-1]  # First n-1 parameters are coefficients for the independent variables
    b = params[-1]  # The last parameter is the intercept
    y_pred = np.dot(x, a) + b
    return np.sum((y_pred - y) ** 2)  # Minimize the sum of squared errors

# Define constraints, for example a1 + a2 + a3 + a4 + a5 = 1
# constraints = ({'type': 'eq', 'fun': lambda params: np.sum(params[:-1]) - 1})
constraints = ()
# Initial guess for parameters
initial_guess = np.random.rand(x.shape[1] + 1)
# Minimize the objective function
result = minimize(objective, initial_guess, args=(x, y), constraints=constraints)
# Get the results of multivariate linear regression with constraints
coefficients = result.x[:-1]  # First n-1 parameters are coefficients for the independent variables
intercept = result.x[-1]  # The last parameter is the intercept

# Output multivariate linear regression results
print(f"Coefficients (Coefficients): {coefficients}")
print(f"Intercept (Intercept): {intercept}")

# Use the model to make predictions
y_pred = np.dot(x, coefficients) + intercept
y_pred1 = np.dot(x1, coefficients) + intercept

# Calculate Mean Squared Error (MSE)
mse = mean_squared_error(y, y_pred)
# Calculate R-squared
r2 = r2_score(y, y_pred)
# Calculate Mean Absolute Error (MAE)
mae = mean_absolute_error(y, y_pred)

# Output model evaluation results
print(f" Mean Squared Error (MSE): {mse}")
print(f"R-squared: {r2}")
print(f" Mean Absolute Error (MAE): {mae}")

# Visualization
years = np.arange(2012, 2012 + len(x))
years1 = np.arange(2012, 2012 + len(x1))
plt.plot(years, y, label='Actual values', color='blue', marker='o')
plt.plot(years1, y_pred1, label='Predicted values', color='red', linestyle='--', marker='x')
plt.xlabel('Year')
plt.ylabel('Pet food market size')
plt.legend()
plt.show()

# 灵敏度分析
variations = [-0.15, -0.10, -0.05,0, 0.05, 0.10, 0.15]
best_a1 = coefficients[0]
plt.figure(figsize=(12, 6))  # 设置图形大小
for var in variations:
    new_a1 = best_a1 * (1 + var)
    new_coefficients = coefficients.copy()
    new_coefficients[0] = new_a1
    y_pred_var = np.dot(x1, new_coefficients) + intercept
    linestyle = '--' if var < 0 else '-'  # 负变化使用虚线，正变化使用实线
    marker = 'o' if var < 0 else 'x'    # 负变化使用圆点，正变化使用叉号
    plt.plot(years1, y_pred_var, label=f'{"+" if var > 0 else ""}{var*100:.0f}%', linestyle=linestyle, marker=marker)

plt.title('Sensitivity Analysis', fontsize=14)  # 修改图表标题
plt.xlabel('Year', fontsize=12)  # 调整x轴标签字体大小
plt.ylabel('Amount of pet', fontsize=12)  # 调整y轴标签字体大小
plt.legend(fontsize=10, title_fontsize=12, loc='upper left')  # 优化图例
plt.grid(True, linestyle='--', alpha=0.5)  # 添加网格线
plt.tight_layout()  # 自动调整子图参数以填充图形区域
plt.show()

# Residual Analysis
residuals = y - y_pred
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.scatter(y_pred, residuals)
plt.axhline(0, color='red', linestyle='--')
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.title('Residuals vs Predicted')

plt.subplot(1, 2, 2)
plt.hist(residuals, bins=20)
plt.xlabel('Residuals')
plt.ylabel('Frequency')
plt.title('Residuals Distribution')
plt.show()

# Use statsmodels to calculate more performance metrics
x_with_intercept = sm.add_constant(x)
model = sm.OLS(y, x_with_intercept).fit()

# Output F value and p-value
f_value = model.fvalue
p_value = model.f_pvalue
print(f"F value: {f_value}")
print(f"p-value: {p_value}")

# Output confidence intervals
conf_intervals = model.conf_int()
print("Confidence intervals:")
print(conf_intervals)

print(y_pred1)

# Output model summary
# print(model.summary())